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Abstract 

■ We investigate the dynamics of two bosons trapped in an infinite one-dimensional optical lattice 

potential within the framework of the Bose-Hubbard model and derive an exact expression for the 
^ i wavefunction at finite time. As initial condition we chose localized atoms that are separated by a 

distance of d lattice sites and carry a center of mass quasi-momentum. An initially localized pair 
[d = 0) is found to be more stable as quantified by the pair probability (probability to find two 
atoms at the same lattice site) when the interaction and/or the center of mass quasi-momentum 
is increased. For initially separated atoms {d ^ 0) there exists an optimal interaction strength 
for pair formation. Simple expressions for the wavefunction, the pair probability and the optimal 
interaction strength for pair formation are computed in the limit of infinite time. Whereas the 
time-dependent wavefunction differs for values of the interaction strength that differ only by the 
sign, important observables like the density and the pair probability do not. With a symmetry 
analysis this behavior is shown to extend to the A^-particle level and to fermionic systems. Our 
results provide a complementary understanding of the recently observed [Winkler et al, Nature 
(London) 441, 853 (2006)] dynamical stability of atom pairs in a repulsively interacting lattice gas. 
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I. INTRODUCTION 



The physics of particles trapped in periodic potentials has been a topic of extensive re- 
search since the early days of quantum theory [l -^. The progress in experimental techniques 
during the last two decades which has led to the realization of atomic Bose-Einstein con- 
densates (BECs) in optical lattice potentials and the possibility to tune the inter-atomic 
interaction via Feshbach resonances has renewed this interest |^-l6|. Because of the high 
control over the system's parameters and the absence of strong dissipation channels, it is 
possible to simulate periodic systems isolated from other effects like, for example, phononic 
degrees of freedom which play an important role in solid state physics. Due to this isola- 
tion not only ground state properties but also the excited states and the dynamics of such 
systems play a crucial role in the understanding of present experiments. 

Of special interest for this work is the experiment recently done by Winkler et al. [3] 
who show that two repulsively interacting atoms initially prepared at one site of an op- 
tical lattice potential separate less rapidly than their non-interacting counterpart. After 
this experiment several theoretical works on the topic followed, most of them within the 
framework of the Bose- Hubbard model. Different two-body problems (with several trapping 
potentials and interactions) are investigated in |8l-ll5|. Other authors study the effect using 



a larger ensemble of particles |16l . Il7j . But pairing induced by a repulsive interaction is 
not restricted to bosonic atoms, it is a relevant topic for fermions as well 18|-|20|. A work 
that does not use the framework of a Hubbard or Bose-Hubbard model is an extension of 
an older work on fermionic pairing in the context of high temperature superconductivity 



by Mahajan and Thyagaraja 



21 



22j, who point out that the effect of pairing by repulsion 



has three ingredients, namely quantum mechanics, a periodic potential and a short range 
interaction. 

In this work we are interested in pairing by repulsion (or attraction) from the viewpoint of 
dynamics. Usually the computation of quantum dynamics is a difficult task and a hot topic 



of actual research 



23 



271]. Because of the complexity of the problem approximations are 



often a must and most computations are based upon numerical methods. Therefore, exactly 
solvable models are of special interest. In this work we present an exact solution of the two- 
body dynamics of two initially localized atoms carrying a center of mass quasi-momentum 
within the framework of the Bose-Hubbard model. The solution is then applied to the 
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physical problem of pairing of two bosonic atoms trapped in an optical lattice potential. 
Common observables are found to depend only on the magnitude of the interaction strength 
but not on its sign. With a symmetry analysis this behavior is shown to extend to the 
corresponding (bosonic or fermionic) A^-particle system. The paper is structured as follows. 
In Sect, iniwe derive the time-dependent wavefunction and discuss results of the dynamics 
at finite time. A simple way to solve the dynamics in the limit of infinite time is presented in 
Sect. Illli The following section (Sect. IIV|) is concerned with the question how the dynamics 
of the pair and some generalization for the A^-particle system depend on the sign of the 
interaction strength. Finally, in Sect. |V]we summarize the findings and give a short outlook. 



II. DYNAMICS AT FINITE TIME 



A. The time-dependent wavefunction 

The Bose- Hubbard model mostly is used in its second quantized version [28(. The Hamil- 
tonian describing particles in an infinite one-dimensional periodic potential reads 

U 



H = -jJ2 {^ik+i + bi^ik) + ^ ^ n„ (n„ - 1) . (1) 



As usual bl^/ba is the creation/annihilation operator and ha the number operator for a par- 
ticle at site a, which is described by a Wannier function centered around the lattice site. 
J and U denote the tunneling rate between neighboring sites and the on-site interaction, 
respectively. When treating two-particle systems within the Bose-Hubbard model its repre- 



sentation in first quantization is more convenient 
equation is then given by the expression 



The two-particle Schrodinger 



- J [A, + A^] ^(x, y) + US.^y^ix, y) = E^{x, y). (2) 

The coordinates x and y are elements of the lattice F = Z (lattice constant a = 1) and for 
the Laplace operator on the lattice we define /S.x'^{x) = \E'(a:; + 1) + \I^(a; — 1). To exploit 
the translational symmetry we introduce relative and center of mass coordinates r = x — y, 
R = ^-T^- The center of mass coordinate is an element of the lattice F' = (TLjl:^^ while the 
relative coordinate is an element of the original lattice F. For the wavefunction the usual 
ansatz \l/(x, y) = e*^'^\E'i^(r) can be made. The time-independent Schrodinger equation thus 



3 



transforms to 



[-Jk^v + U5r,^] ^x(r) = EK-^Kir) for K e [-tt, tt] , (3) 



where Jk = 2J cos{K/2) is an effective hopping parameter. Please see [4l|] for an argument 
why K G [— TT, tt] although the center of mass lattice T' has a lattice constant of acM = |- 
The solutions of the above eigenvalue equation consist of one bound state and a continuum of 
scattering states for each value oi K 7|,|8|]. In contrast to the corresponding continuum model 
with delta function interaction [si!] , the bound state exists also for repulsive interaction. For 
the bound state wavefunction and energy there are two distinct expressions, one for attractive 
interaction and one for repulsive interaction 
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(f/i^-yZ4TT)''',for t/>0 



, Vf7^^4J|, for f/ < 0, 
Ebs{K) = { I ^' (5) 

+ 4J|, for f/ > 0. 

Uk = U / {2Jx) denotes an effective interaction parameter. We note that the bound state 
has the same density p^^{R,r) = |\E'^'^(i?, r)|^ for repulsive and for attractive interaction. 
Its energy lies below the scattering continuum for U < (lowest state in energy) and above 
it for U > (highest state in energy). The scattering states appear in the literature in an 
unnormalized form. We provide here the normalized functions 



^UkAR, r) = , ^ <^ cos(A;r) + —f- sm{k \r\) \ , (6) 



sin(A;)^ 

E{K,k) = -2JKCOs{k), (7) 

k G [— TT, 7r]. Note the symmetry of the spectrum for K — )■ —K and k — i- —k. 

Given an initial state at t = the dynamics of the system can be computed with the 
time evolution operator U{t) = e~*^* (^ = !)• Making an expansion in the eigenbasis of 
the Hamiltonian this calculation can be done explicitly. As an initial condition we make 
the choice ^ofl{R,r) = ^^^r,o to describe two particles sitting at the same lattice site 
and having a center of mass quasi-momentum Q, or \E'o,d(-R, r) = ^^^^ (6r,d + 6r-d) for 
two particles that are initially separated by a distance of d lattice sites. Expanding the 



initial state into the above basis and propagating it in time we find the time-dependent 
wavefunction 

piQR 



^u<o{R,r,t) 



'2ti 
2^ 



Cdl{r, t) + Cd 



Cdl{r, t) + Cd 



\Uc 



\r\ + \d\ 



\r\ + \d\ 



(8) 



(9) 



Cd is a normahzation constant which equals one for c? = and \/2 for d 7^ 0. The function 
/(r, t) is given by 

dk 



I{r,t) 



' ^i2jQt cos 



vr 



fn{k) = cos{kn) + 



Ur 



sin(/c) 



ik) fd{k)fr{k) 

1 + 

sm{k \n\), 



(10) 



(11) 



where n E'L. It is a closed form integral expression that can easily be evaluated numerically 
with standard computer algebra programs. In what follows we use the notation 



vl>(i?,r,t) 



(12) 



to have an explicit expression for the time-dependent wavefunction in the relative coordinate 
at hand. We note that <l>Q(r, t) stands for the expression in brackets in Eq. ^ [Eq. ([9])] 
for U < (t/ > 0). Since the wavefunction in the center of mass coordinate is a spa- 
tially oscillating exponential function, the two-particle density depends only on the relative 
coordinate. 



B. Dynamics of the pair probability and the density in the relative coordinate 

We are interested in the question how a pair forms or dissociates in time. To quantify 
the process we use the pair probability which measures the probability to find two particles 
at the same lattice site 15|. The pair probability is of special interest because it is one 
of the observables used in the experiment on repulsively bound atom pairs [7]. For a two- 
particle system it can be defined as the expectation value of the pair operator Ppair = 
Ylia& l^a — 2) ('^a = 2| which in first quantization reads Ppair = ^rfl- Therefore, the pair 
probability equals the density in the relative coordinate at point r = 0, 

Ppa^r{t)=\<^Q{r = Q,t)\\ (13) 



Our second observable is the variance of the pair probabihty which we use to measure 
its fluctuations. We note that the pair operator is a projector and fulfills the relation 
Ppair = Ppair- Hencc, the variance of the pair probability can be written as a function of the 
pair probability itself: 

^pairi^) ~ (^(^Ppair ~ (^Ppair^^ ^ = Ppairif) [l — Ppairif)\ ■ (14) 

As a third observable we monitor the density in the relative coordinate in order to get an 
impression how the initial wave packet changes its shape. It is given by 

p(r,t) = |<|.Q(r,^)|^ (15) 

All calculations have been carried out for J = 1 and Q = 0. Other values of the hopping 
parameter and the center of mass quasi-momentum scale the time t linearly with Jx = 
2J cos{K/2) and the on-site interaction U inversely with 

The three computed quantities, Ppairit), crlairi't) p(r, t) depend only on the magnitude 
of the interaction strength and not on its sign. Therefore, it is sufficient to treat the case 
U > 0. We note that this behavior is a nice explanation for the existence of repulsively 
bound atom pairs. Using the initial conditions from Sect. [Tll the repulsively interacting 
particles act as if they would attract each other and vice versa. Hence, there is a strong 
relation between the phenomena of binding by repulsion and binding by attraction. In 
Sect. II VI we will derive a similar result for the A^-particle system which gives a direct link to 
the experiment on repulsively bound atom pairs. 

When we start the dynamics with two atoms localized at the same lattice site [see 
Fig. l(a)| the pair probability starts at Ppairit = 0) = 1. For f/ = it oscillates for small 



times and then rapidly goes to zero as 1/t for large t. For U ^ this behavior changes. 
Again, we have some oscillations for small times but for large t the pair probability evolves 
towards a finite asymptotic value. This asymptotic value is reached more quickly in case of 
a larger interaction strength. Additionally, the asymptotic value of the pair probability is 
found to increase with the interaction strength. For an arbitrary fixed time the pair prob- 
ability is an increasing function of the interaction strength. Hence, we conclude that if we 
start the dynamics with a pair (two particles at the same lattice site), a large interaction 
strength stabilizes the pair. Since the effective interaction strength scales inversely with 
cos((5/2) [42] also a large center of mass quasi-momentum stabilizes an initially prepared 
atom pair (for U 0). 



If we start the dynamics with two atoms sitting at adjacent lattice sites [d = 1, Ppair{t = 
0) = 0] the quahtative shape of the pair probabihty changes although some phenomenology 
stays the same [see Fig. [T(b)] . In the non-interacting case {U = 0) we again have some 
oscillations for small times, and again Ppair{t) goes to zero for large times. Additionally, like 
for d = the pair probability evolves towards a finite asymptotic value as long as U 7^ 0. 
An important difference is that it has a maximum as a function of U for all finite t (except 
for very small times). The existence of this maximum can be understood as an interplay 
between two effects in the process of pair creation. As long as the two particles are separated, 
a small interaction strength is needed for pair formation. But when the particles have come 
together and are sitting at the same lattice site, a large interaction strength is needed to 
stabilize the pair. Hence, there exists an optimal interaction strength for pair formation 
which we have computed analytically in the limit of infinite time, see Sect. Illli 

We note that the argumentation used to explain the existence of an optimal interaction 
strength for pair formation holds for repulsive as well as for attractive interaction. In the 
latter case a strong attractive interaction keeps the atoms apart from each other which is 
counterintuitive at first sight. The effect can be understood in the following way. Two 
atoms coming together at one lattice site lower their interaction energy and due to energy 
conservation they gain additional kinetic energy. When the width of the first Bloch band 
(4Jq) is smaller compared to U this process is strongly suppressed. 

The variance of the pair probability as a function of the pair probability itself CTpairit) = 
Ppair{t) [1 — Ppairit)] IS Symmetric around Ppair{t) = \ where it takes its maximum; and for 
Ppair{t) = 0, 1 it equals zero. We conclude that the variance is small whenever the pair 
probability is near to zero or one. Large fluctuations occur only in case of intermediate 
pair probabilities. For d = and d = 1 the variance of the pair probability plotted as a 
function of U and t has very similar characteristics, see Fig. [2l For U = both curves start 
at CTpairit) = 0, then they quickly increase, oscillate and go to zero for large times. Like in 
the case of the pair probability, this behavior changes when a finite interaction strength is 
turned on. The two surfaces now converge to an asymptotic value larger than zero. For all 
fixed and not too small times t (t > 2 is sufficient), both functions possess a local maximum. 
In case of d = 1 this maximum coincides with the maximum of the pair probability since 
the latter mentioned is smaller than ^ for all U and t, and hence in this regime a large pair 
probability always goes hand in hand with large fluctuations of this quantity. 
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What about the density? If we start the dynamics with two particles at one lattice site 
[d = 0) and set U = 5 the initially localized wave packet splits into three wave packets 
that are propagating in time, see Fig. [3J The wave packet centered around r = describes 
the time evolution of a pair while the other two wave packets describe the separation of 
the two particles. We remind that the density at r = equals the pair probability. The 
fact that we see two wave packets describing the separation of the particles is due to the 
bosonic symmetry of the wavefunction. For two particles that initially sit at adjacent lattice 
sites {d = 1) the picture qualitatively looks the same, see Fig. HI As one would expect the 
wave packet describing the time evolution of the pair is less pronounced which means that 
the particles separate with a higher probability. If we start the dynamics with larger initial 
distances the dynamics of the density becomes more complicated. For d = 5 the initial wave 
packet splits into four wave packages that have a much broader shape than in the case of 
d = or d = 1, see Fig. [51 The wave packet describing the pair propagation has nearly 
vanished. 



III. ASYMPTOTIC DYNAMICS 



In the previous section we have seen that the pair probability is converging rapidly towards 
an asymptotic value for large times. We are interested in this asymptotic value and therefore 
compute the wavefunction, the pair probability and its variance in the limit of infinite time. 
The calculation can easily be done with the help of the Riemann-Lebesgue lemma which is 
used to show that the function I{r,t) defined in Eq. ( ITOj) vanishes in the desired limit, see 
App. S The asymptotic wavefunction reads 

iQR \ xlr-l+MI 

n<o{R.r) = c, ^ 1^^ (Ji + UI-\Uq\) , (16) 



27r . /i I rr2 



i + ui 

n>o{R. r) = (Uq - Jul + 1 ) . (17) 



2^^M + 1 



We note that it is not normalized anymore. This is because \l/(t) converges to only 
pointwise and not in the P norm. The same behavior can be found for example for a 



textbook Gaussian wavefunction describing the dynamics of a particle in free space 32|. For 
all times the wavefunction is normalized. In the limit of infinite time it converges pointwise 
to zero but it does not converge to the constant zero function in the norm. In our case this 
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property can make the calculation of other asymptotic quantities more demanding because 
one may not be allowed to interchange summations coming from the scalar product with 
the time limit. 

/^From the asymptotic wavefunction we compute the asymptotic pair probability 

U}, / I \2MI 

Pvair = cljJ^[^l + Ul-\UQ\) , (18) 

where Un = tM- = ti — tttavt- For d = the formula reduces to Pil„v = which is 

^ 2Jq 4Jcos(Q/2) pa.tr Uq+1 

a strictly increasing function of the effective interaction parameter and approaches one for 



large effective interaction strengths, see Fig. 6(a) Hence, also from the asymptotic pair 
probability we can conclude that a large interaction strength or a large center of mass quasi- 
momentum stabilizes an initially prepared pair. As expected, the asymptotic pair probability 



has a local maximum for all ci 7^ 0, see Fig. 6(a) for d = 1, 2. The location of this maximum 
defines the optimal interaction strength for pair formation which is given by 



T jmax I 



1 + ^ 



\d) = ±^ ^ . (19) 

We note that UQ°'^{d) is a strictly decreasing function of the initial distance d. Hence, if we 
chose a larger initial distance the optimal effective interaction strength for pair formation 
will be smaller. 

As we have seen in Sect. [TTl the variance of the pair probability is a function of the pair 
probability itself. The same holds for the asymptotic variance of the pair probability which 

reads 

at.,. = P;,,, [1 - P;,„] , (20) 



pair 



see Eq. (fT^ . It has a local maximum for all values of the initial distance d, see Fig. 6(b) 



for d = 0,1,2. In case of > the asymptotic pair probability is always smaller than ^ 
and we can conclude that the maxima of the asymptotic pair probability and its variance 
coincide. Hence, in this regime a large pair probability always leads to a large variance of 
this quantity. 



IV. DYNAMICAL SYMMETRY IN THE iV-PARTICLE SYSTEM 



In Sect. Ill Bl we have shown that the time-dependent pair probability, its variance and 
the time-dependent density do not depend on the sign of the interaction strength U. This 
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behavior is an important finding because it provides an explanation for tlie existence of 



repulsively bound atom pairs that is complementary to the one given in p']. If a repulsive 
interaction leads to the same time-dependent pair probability as an attractive interaction 
of the same magnitude, why should there be no dynamical stability of an initially prepared 
pair? Nevertheless, until now we have investigated the dynamics of two particles whereas 
the experiment has been performed with a gas consisting of approximately 2 ■ 10^ atoms. 
To be able to make statements also in this regime, we extend the result on the invariance 
of the three observables under a change of the sign of U from Sect. IIIBI and show that a 
similar statement holds for the A^-particle system, provided one chooses the "right" initial 
conditions. Interestingly, for its proof one does not need to specify the statistics of the 
particles - it works for bosons and for fermions alike. As a prerequisite we discuss simple 
relations between the spectra and eigenf unctions of the attractive and the repulsive Bose- 
or Fermi-Hubbard model. 

It is well known that the Bose-Hubbard model with two lattice sites (Bose-Hubbard 
dimer) possesses a symmetry connecting the attractive and the repulsive Hamiltonian. This 
symmetry leads to relations between static proper ties like for example the energy spectrum 



of the attractive and the repulsive system 



33N36|. As was first realized and quantified in 



37l | the symmetry also affects the dynamics and leads to a dynamical symmetry in the 
Bose-Hubbard dimer. The authors could show that the Bose-Hubbard model dictates an 
equivalence between the time evolution of the survival probability and fragmentation of 
the attractive and the repulsive system if all particles initially reside in one of the two 
wells. A short time afterwards it was shown that the time evolution of expectation values 
in the Fermi-Hubbard model under certain conditions does not depend on the sign of the 



interaction strength 
some of the results to be found in 
conditions. 



The main result we present here (Theorem 2) is an extension of 
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for a more general class of operators and initial 



A. A relation between the spectra and eigenfunctions of the attractive and the 
repulsive model 



We start the analysis with the definition of the unitary operator [3j 

R={a^^{-ir a^}. (21) 
10 



It changes the sign of the creation and the annihilation operator at every second lattice site 
and thereby shifts the quasi-momentum of each particle by tt, see App. [Bl By (ia we denote 
a bosonic or fermionic creation operator for a particle at lattice site a. Possible spin indexes 
are surpressed since they do not play a role. For convenience we assume a one-dimensional 
infinite lattice F = Z but everything still works for a finite lattice with periodic boundary 
conditions for an even number of lattice sites and in higher dimensions. When we assume 
the Bose-Hubbard or Fermi- Hubbard Hamiltonian j^sl, the following relation holds js^] 

RH{U)R = -H{-U). (22) 

Using this equality it is easy to proof a statement about the spectrum and the eigenfunctions 
of the Hamiltonian H{U). 

Theorem 1: Let |\E'(f/)) be an eigenf unction of the Hamiltonian H{U) with eigenvalue 
E{U). Then \^{—U)) = R\'if{U)) is an eigenfunctions ofH{—U) with eigenvalue E{—U) = 
-E{U). 

Proof: Assume we are given the solution of the time-independent Schrodinger equation 
with interaction strength t/, H{U) |^(f/)) = E{U) \'^{U)). When we let R act on both sides 
of this equation we find H{-U)R \^{U)) = -E{U)R \^{U)). ■ 

Hence, the simple relation between the Hamiltonian of the attractive and the repulsive 
model leads as well to simple relations between their spectra and eigenfunctions. Theo- 
rem 1 can also be interpreted as a statement about the existence of repulsively bound states 
in the two models which can be seen as follows. We know that the wavefunctions of the 
attractive and the repulsive system are related by the operator R. Translated to first quan- 
tization this relation reads \E'(a;i, xat, — f/) = {—1Y^^'''^^^'^{xi,...,xn-,U) (see App. [Bl). 
and consequently the density of '^{U) and of '^{—U) equal each other. From this it fol- 
lows that if \l/(f/) is square summable the same must be true for \l/(— f/). Therefore, the 
Bose-Hubbard and the Fermi-Hubbard model with attractive and repulsive interaction of 
the same magnitude have an equal number of bound states (square summable wavefunc- 
tion). This surprising finding is a pure lattice effect and contrasts with quantum mechanics 
in a continuous coordinate space where potentials usually do not have bound states anymore 
when their character is changed from attractive to repulsive. 

As an example we mention the solution of the time-independent Schrodinger equation 
[Eq. ([3])] from Sect. Ill A[ As already mentioned the spectrum of the Hamiltonian consists 
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of a scattering continuum [Eq. ([7])] whose energy is invariant under the operation E — )■ —E 
and a bound state below or above the scattering continuum for attractive and repulsive 
interaction, respectively [Eq. (^]. The bound state wavefunctions of the two systems differ 
by a factor of (-l)!*^! = (-1)^+*'. 

B. Invariance of time-dependent expectation values under the transforma- 
tion U ^ -U 

As we have seen in the previous subsection, Eq. fl2T]) leads to simple relations between the 
spectra and eigenfunctions of the attractive and the repulsive system. This certainly affects 
also the dynamics. Upon operation of R, the transformation of the spectrum {E — )■ —E) 
leads to a reversal of time; the transformation of the wavefunction (\E' — )■ R"^) is not that 
intuitive, see App. [Bl Nevertheless, with a little more effort needed to derive the time- 
independent result it is possible to make a statement about time-dependent properties as 
well. We start with a definition. 

Definition: We call an operator O real if it has only real coefficients when being expressed 
with the creation and annihilation operators da and d]^. An operator that is given by matrix 
elements in the Wannier basis [44] called real if all its matrix elements are real. 
Using the definition we state the main result of this section: 

Theorem 2: Let O be a hermitian operator that can be written as O = Or + iOi, where Or 
and Oi are real operators and fulfill the relations R Or R = Or and R Oi R = —Oi. We 
assume that H{U) is the Bose-Hubbard or Fermi-Hubbard Hamiltonian. The initial condition 
reads = if + ix where (p and x are assumed to be real functions on the lattice = Z^. 
Additionally, (p and x are eigenfunctions of R with different eigenvalues, in formulas R\(p) = 
± Iv?) and R\x) = -F|x)- Then the following relation for the time- dependent expectation 
value of 6 IS true: 0{U,t) = 0{-U,t). 

Proof: We give here the proof only for the special case when \E'o is a real function on 
and O is a real operator in analogy to The general proof can be found in App. O We 
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write the expectation value of O as 

0(U,t) = (^ol e'^^'O e-'^^' |^o) (23) 
= (^o|cos(^c/t)dcos(iJ[/t) l^o) 
+ (^ol sm{Hut)d sm{Hut) |^o) 
- 2 Im (*o| sm{Hut)d cos{Hut) |^o) 

The problem we have to manage is to change the sign of U with an insertion of = 1 
factors without changing the direction of time. This is because Re~^^^^R = e*^~^*. On 
that account, let us analyze the matrix element (\E'o| sin{Hut)0 cos{Hut) {"^o). All operators 
have only real coefficients when being expressed with cIq, and aj^, or when being expanded 
in the Wannier basis. Additionally, the creation and annihilation operators produce only 
real numbers when acting on occupation number states (or Wannier basis states which is 
the same in this setting) and the wave function at time zero \E'o can be expanded into the 
Wannier basis with real coefficients only (it is a real function on F^). Hence, the overall 
matrix element is real and its contribution to the expectation value of O vanishes. We find 

0{U, t) = (^ol cos{Hut)d cos{Hut) + sm{Hut)d sm{Hut) |^o) • (24) 

Now we can insert R"^ = 1 factors between all operators and thereby change the sign of the 
interaction strength U without changing the direction of time. When we additionally use 
the two relations ROR = O and R'^o = i^E'o the result of Theorem 2 can be shown. ■ 

Let us discuss the applications of Theorem 2. First we note that it applies to the two- 
particle dynamics of Sect. [TTl The initial condition '^o,o{R,''^) = ^^5r,o reads l^&o.o) = 



Ex,yer ^'^^'^' ^^.^^t&t |o) = J2^^r (^i) 1^) second quantization, also see App. M 
When we let R act on this state we find -R|\E'o) = l^o)? but obviously the expansion co- 
efficients are not real. This can be circumvented by using the Hamiltonian of Eq. ([3]) as 
the starting point. In second quantization it reads H{U) = —Jr J2aer ^Ifia+i + h.c. + Utiq 
and fulfills Eq. (!22|) as well. Restricted to the one-particle subspace (of the Fock space) 
its dynamics are the one governed by the original Hamiltonian in the relative coordinate 
(the wavefunction in the center of mass coordinate is constant in time), and hence as initial 
condition we have to use \l/o Q(r) = 6r,o- In second quantization this reads |^E'o,o) = 
still fulfills the relation R |\E'o) = l^o); and additionally has real expansion coefficients. The 
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same calculation can be done with the initial condition "^Q^^ir) = ^ {Sr^ + Sr-d) which in 
second quantization reads I'^od) = {h\ + |0). Acting with R on the state we find 
R\^o,d) = (—1)'^ |^o,d)- Hence, \l/og and \E'q^ qualify as possible initial conditions. As ob- 
servables we used the time-dependent pair probability, its variance and the time-dependent 
density. For their computation we have to evaluate expectation values of the operators 
Ppair = Ylaer I'^a = 2) = 2| which, whcu acting on one-particle states in the relative 
coordinate, can be written as Ppair = ^r=o^r=o and fij. = h\.hr. Obviously, Ppair and fir 
both are real and fulfill the relations RPpairR = Ppair and RfirR = fir. Therefore, all three 
observables and with a little trick also the initial conditions qualify for an application of 
Theorem 2. The result explains why the dynamics of our observables does not depend on 
the sign of the interaction strength U . 

How can Theorem 2 be used to extend the invariance properties of the pair probability, 
its variance and of the density to the A^-particle system? First, we need to find appropriate 
observables. The pair probability (probability to find two particles at the same lattice site) 
br N particles can be written as the expectation value of a pair operator as well. It reads 

2 . 

Ppair = ^y^,\'^n) mn (v?n| , (25) 
neN 

where by {v^nj^g^ we denote the Wannier basis of the A^-particle Hilbert space [44]. The 
number m„ counts the pairs in each Wannier basis state and the prefactor ^ assures the nor- 
malization. We highlight that the pair operator is a A^-particle operator, and consequently 
the full wavefunction is needed to compute its expectation value. The fact that the Wannier 



basis states are eigenf unctions of ^ 45] can be used to show that R Ppair R = Ppair holds. 



For the computation of the variance of the pair probability we also need to compute the 
expectation value of Ppair which has the same properties. As an equivalent for the density 
in the two-particle system we choose the p-particle density (p < A^) that is given by the 
expectation value of the operator 

fiP{xi,...,Xp) = al^ ■ ... -ai^a^p ■ ... ■ a^^- (26) 

Obviously, it qualifies for an application of Theorem 2 as well. As initial states for the two- 
particle dynamics we have chosen particles that (in the relative coordinate) are described by 
Wannier basis states, and on account of this describe particles that are localized to single 
lattice sites. Since Wannier basis states describe real functions on and are eigenfunctions 
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of R [45], we can choose them as possible initial conditions for the A'"-particle dynamics as 
well. Nevertheless, Theorem 2 tells us that the class of all possible initial states is much 
larger. 

We conclude that Theorem 2 generalizes the invariance properties of the pair probability, 
its variance and the density under the transformation U — )■ —U to a system with N bosonic 
or fermionic particles. Additionally, Theorem 1 tells us that repulsively bound states are 
not a speciality of the two-particle system but will always occur when the attractive system 
has bound states. The experiment on repulsively bound atom pairs has been performed 
with atom pairs that initially have been localized to single lattice sites. Whether these 
initial states qualify for an application of Theorem 2 is difficult to say. Nevertheless, if we 
approximate the experimental initial state with a pure Wannier basis state we can show 
that the dynamics do not depend on the sign of the interaction strength U . The stability of 
repulsively interacting pairs therefore becomes very intuitive when we understand that the 
lattice structure of the coordinate space forces the attractive and the repulsive system to 
act very similar. 

V. SUMMARY AND OUTLOOK 

To summarize the findings, we have computed an exact expression for the time-dependent 
wavefunction for two bosons trapped in an infinite one-dimensional optical lattice potential 
within the framework of the Bose-Hubbard model. As initial conditions we have chosen 
localized atoms that are separated by a distance of d lattice sites and carry a center of mass 
quasi-momentum. An initially localized pair (d = 0) is found to be more stable as quantified 
by the pair probability when the interaction and/or the center of mass quasi-momentum is 
increased. In contrast, for two initially separated atoms there exists an optimal interaction 
strength for pair formation. 

To gain further information we have monitored the variance of the pair probability and the 
density in the relative coordinate during the dynamical process. Analytical expressions for 
the wavefunction, the pair probability and the optimal interaction strength for pair formation 
have been derived in the limit of infinite time. We had to give two distinct expressions for the 
time-dependent wavefunction for positive and negative interaction strength. In contrast, the 
pair probability, its variance and the density in the relative coordinate are invariant under 
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the transformation U ^ —U. This leads to the conclusion that for our initial conditions 
the three observables have the same dynamics when being propagated with the attractive 
or with the repulsive Hamiltonian. 

In the second part of the paper we have extended this result and shown that also in the A^- 
particle system there exist time- dependent observables that stay the same when the sign of 
the interaction strength U is changed. The time-dependent pair probability of the iV-particle 
system, its variance and the time- dependent p-particle density {p < N) belong to this class. 
Additionally, we have discussed a simple relation between the spectra and eigenfunctions 
of the attractive and the repulsive Bose- or Fermi-Hubbard model. By showing that the 
dynamics of the pair probability is the same in the attractive and in the repulsive A^-particle 
system, we provide a complementary understanding for the recently observed ^] dynamical 
stability of atom pairs in a repulsively interacting lattice gas. 

The explicit expression for the time-dependent wavefunction we have computed allows 
one for studies which go far beyond the scope of this work. It would be an interesting 
and challenging task to compute for example the one-particle reduced density matrix and 
to study questions of entanglement during the dynamical process. Additionally, there is 
evidence that the computation of the exact time-dependent wavefunction on the infinite 
Bose-Hubbard lattice is possible for other classes of initial conditions, too. 

The finding that the Bose- or Fermi-Hubbard Hamiltonian with attractive and repulsive 
interaction have an equal number of bound states suggests that there is hope to experi- 
mentally find also repulsively bound states consisting of three jsoj or possibly even more 
particles. In addition, the result on the equivalence of expectation values under the trans- 
formation U — )■ —U can be applied to the probability to find m {m < N) particles at one 
lattice site as well. The same holds for the probability to find two particles at one lattice 
site and another one at a neighboring lattice site, which would model the situation of a 
bound state consisting of a dimer and a monomer. The latter suggests that even though 
such objects lie energetically above the scattering continuum, they nevertheless should be 
dynamically stable when being part of a dilute lattice gas. 

Financial support by the DFG is gratefully acknowledged. 
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Appendix A: Asymptotic behavior of I(r, t) 



To compute the asymptotic behavior of the wave function, the pair probabihty and its 
variance, we need to compute the infinite time hmit of the function /(r, t), see Eq. f llUp . This 
can be done with the help the Riemann-Lebeserue Lemma which we state in the following 
form [40;]: 

Riemann-Lebesgue Lemma: Let / be Lebesgue-integrable on [— 7r,7r], then 

fix) e'^'^dx = 0. (Al) 

-TT 

In order to bring /(r, t) to a form that the Lemma is applicable, we do the coordinate 



transformation q = cos{k). Using the relation sin[arccos(g)] = y 1 — we find 

dq e^^'^o*'? 

Hr,t) = ^-5— /rf[arccos(g)]/4arccos(g)] (A2) 

with = cos(fcra) + sin(fc |?t,|). We note that the function /„[arccos(g)] is continuous 

for all n E Z. The expression ^ — ^ — in the denominator of Eq. (IA2p is continuous as 

well, and hence the Riemann-Lebesgue Lemma is applicable, leading to the result 

limJ(r,t) = 0. (A3) 



t—>-oo 



Therefore, the contribution from the scattering states to the pair probability and to the 
density vanishes in the limit of infinite time. 



Appendix B: The action of the operator R 

In this appendix we provide formulas for the action of the operator R defined in Eq. (!2T|) 
on the annihilation operator of a particle with quasi-momentum k and on the wavefunction 
in first quantization in coordinate and in quasi-momentum space. Since the wavefunction in 
second quantization can be expressed with creation and annihilation operators, the action 
of R on it is obvious. To compute its action on the wavefunction in first quantization we 
recall the relation between first and second quantized wavefunctions 

|vi>)= ^{xu...,x^)dl...dijo). (Bl) 

3;i,...,a:jver 
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Here denotes the wavefunction in second quantization, xa?) the wavefunc- 

tion in first quantization and |0) is the vacuum state. As aheady defined in Sect. IIVI 
the operator dj. denotes a bosonic or fermionic creation operator for a particle at site 
X G r. Possible spin indexes are suppressed. When acting with R on |\E') we find 

^1^) = Exi,...,xiver^(^i' •••'^Af)(-1)'''^'"'^'''^ ^^'^ ^Q^ce the action of R on 

the first quantized wavefunction is given by 

x^v) = x^). (B2) 

As one could expect R changes its sign at every second lattice site. For two (or more) 
particles this happens in a chess-pattern-like way. 

How does R act on functions in momentum space? The annihilation operator of a particle 
with quasi-momentum k reads d^ = -^^x&^^'^^^x- Using the relation RdxR = { — lYdx 
we find 

RdkR = cik+n- (B3) 

Having Eq. f lB3|) at hand we can compute the action of R on the first quantized wavefunction 
in quasi-momentum space '^{ki, k^). When we write the second quantized wavefunction 
as |\E') = "^^=iv J^^dki... J^^dkN'^{ki, k]^) d|,^...d|,^ |0) and act with _R on it we find 

M{ki, k^) = ^{ki + TV, kN + 7r). (B4) 

Hence, the operator R shifts the quasi-momentum of each particle by an amount of +7r. 



Appendix C: Proof of Theorem 2 

In Sect. IIVBI we have given the proof of Theorem 2 only for the special case when \E'o 
is a real function on and O is a real operator. Here we assume the general scenario of 
Theorem 2. The initial condition reads = f + '^X^ where ip and x are real functions on F^ 
and eigenfunctions of R with different eigenvalues, in formulas R\(p) = ± Iv'), R\x) = -F 
The observable reads O = Or + iOi with real operators Or and Oj. Additionally, we assume 
that they fulfill the relations R Or R = Or and R Oi R = —Oi. The expectation value of 

d{t) = e'-^vtQQ-iHut ^gg^^g 

0(f/, t) = {ifl d{t) 1^) + (xl d{t) \x) - 2 Im [(^1 Oit) \x)\ . (CI) 
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Let us have a look at the first term of Eq. fICip . It reads {ip\0{t) \ip) = {ip\Or(t) \ip) + 
{ip \ iOi{t) \(p). The term {{p \ Or(t) \lp) is invariant under a change of the sign of the interaction 
strength U because (y9 is a real function and Or is a real operator, see Sect. IIVBI The other 
term can be written as 



{ip\ iOi{t) = i {lp\ cos{Hut)Oi cos{Hut) + sm{Hut)OiSm{Hut) \(p) 



(C2) 



-2 Re 



(v^l sm{Hut)di cos{Hut) \ip) 



Since Oj is a real operator the first matrix element is real. This is because all operators have 
only real coefficients when being expressed with and dj^ or when being expanded in the 
Wannier basis. The creation and annihilation operators produce only real numbers when 
acting on occupation number states (or Wannier basis states which is the same in this setting) 
and \(p) can be expanded into the Wannier basis with real coefficients only (it is a real function 
on r^). The matrix element is multiplied by the imaginary unit i, and hence its contribution 
to the expectation value 0{U,t) is purely imaginary. Because 0{U,t) is a real function this 
contribution has to vanish, thus {(f\cos{Hut)OiCos{Hut) + sm{Hut)OiSm{Hut) {(f) = 0. 
The matrix element in the last term of Eq. ( ]C2p is real because of the same reasons which 
leads to 

(^1 idi{t) \ip) = -2 {ip\ smiHut)di cos{Hut) \ip) . (C3) 
Now we can again insert /^^ = 1 factors between all operators and find 



{^\ idi{t) \^) = -2 ((^1 Rsm{H.ut)RdiRcos{H.ut)R |<^) 



(C4) 



/^From Rsm{Hut)R = — sm{H^ut) we get a factor of (—1) which cancels with the one 
coming from ROiR = —Oi. Similar analysis holds for the second term in Eq. fICip . We 
conclude that {if\ 0(t) \lp) and (x| 0(t) \x) have the wanted invariance property. 
What remains to check is the last term of Eq. ( IClj) . It reads 



Im 



(^1 0{t) \x) 



Im 



+ Im 



{v\^0,{t)\x) 



The first term of the right hand side of this equation can be written as 



Im 



{V\0r{t)\x) 



Im {ip\ cos{Hut)Or cos{Hut) + sm{Hut)Or sin{Hut) \x) 
+i {ip\ sin{Hut)Or cos{Hut) — cos{Hut)Or sin{Hut) \x) 
{(p\ sin{Hut)Or cos{Hut) — cos{Hut)Or sm{Hijt) \x) ■ 



(C5) 



(C6) 



19 



To derive the result we have used the same argumentation as above to show that the matrix 
elements are real. We insert the = 1 factors to change the sign of U and thereby produce 
a factor of (—1) coming from Rsin{Hut)R = — sm{H^ut) which cancels with another one 
coming either from R\ip) = ± |v^) or from R\x) = T Using the same arguments, the 
second term of Eq. (ICSp can be written as 



Im 



{ip\ idi{t) \x) = {^\ cos{Hut)d, cosiHut) + sm{Hut)di sm{Hut) \x) ■ 



(C7) 



Insertion of the obligatory R^ = 1 factors concludes the proof. 
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strength 10 Strength 



(a) (b) 

Figure 1: (Color online) Pair probability {Ppair) as a function of the interaction strength (U) and 
time (t) for J = 1, K = 0, d = (a) and d = 1 (b). With increasing time Ppair is converging 
rapidly towards an asymptotic value. For d = the pair probability increases with the interaction 
strength while for d = 1 there exists an optimal interaction strength for pair formation (except for 
very short times). See Sect. HTBl for more details. All quantities are dimensionless. 
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(a) (b) 

Figure 2: (Color online) Variance {(Tpair) of the pair probability as a function of the interaction 
strength ([/) and time (t) for J = 1, K = 0, d = (a) and d = 1 (b). Both curves have very similar 
characteristics, namely a local maximum for intermediate interaction strengths and the convergence 
towards a constant value for large times. For d = 1 the maximum of the pair probability and its 
variance coincide, see Sect. Ill B] for more details. Therefore, in this regime a large value of the pair 
probability goes hand in hand with large fluctuations. All quantities are dimensionless. 
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Time 




Relative Coordinate 



Figure 3: (Color online) Density [p(r, t)] in the relative coordinate as a function of time (t) for 
U = 5, J = 1, K = and d = 0. It can be seen that there are three wave-packets propagating in 
time. The one propagating to the left and the one propagating to the right (with equal amplitude 
due to the bosonic symmetry of the wavefunction) describe the separation of the two particles. 
The one with the largest amplitude is centered around the origin and describes the time-evolution 
of the pair. See Sect. Ill Bl for more details. All quantities are dimensionless. 
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Relative Coordinate 



Figure 4: (Color online) Same as Fig. [3] but for d = 1. The wave-packet in the middle is less 
pronounced than for d = 0. Hence, there is a smaller probability to find a pair. See Sect. Ill Bl for 
more details. All quantities are dimensionless. 
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Relative Coordinate 



Figure 5: (Color online) Same as Fig. [3] but for d = 5. In the dynamics with a larger initial distance 
new effects like a second splitting of the wave-packet appear. See Sect. IIIBI for more details. All 
quantities are dimensionless. 
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(a) (b) 

Figure 6: (Color online) Asymptotic pair probability [-Ppaj,,] (^) ™d variance [Cp^j^] (b) as a 
function of the effective interaction strength Uq for J = 1 and c? = [dashed blue (highest) curve], 
d = 1 [solid violet (middle) curve], d = 2 [dot-dashed brown (lowest) curve]. The asymptotic pair 
probability has a local maximum for d 7^ 0, the asymptotic variance for all initial distances. Since 
for d 7^ the asymptotic pair probability is always smaller than ^ its maximum and the maximum 
of its variance coincide, see Sect. IIIII for more details. All quantities are dimensionless. 
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